# export LD_LIBRARY_PATH=/usr/lib/jvm/default-java/jre/lib/amd64:/usr/lib/jvm/default-java/jre/lib/amd64/server

NETLOGO_PATH <- '/home/rick2/Apps/netlogo-5.1.0'
MODEL_PATH <- 'hivsim0.nlogo'
NSim <- 10
NTicks <- 1000  # FIXME: tem que ser repetido dentro de simfun()

tic <- Sys.time()

library(parallel)
processors <- detectCores()
cl <- makeCluster(processors)


# initialization
presim <- function(dummy, gui, nl.path, model.path) {
	library(RNetLogo)
	NLStart(nl.path, gui=gui)
	NLLoadModel(model.path)
}

# the simulation
simfun <- function(x) {
	print("running simulation")
	NTicks <- 1000  # FIXME: pq?!
	NLCommand('setup')
	d <- data.frame(hiv=numeric(NTicks), ab=numeric(NTicks), th=numeric(NTicks))
	for(t in 1:NTicks) {
		hiv <- NLReport('count hiv')
		ab  <- NLReport('count ab')
		th <- NLReport('count thCell + count memTh + count effTh')
		d[t,] <- c(hiv,ab,th)
		NLCommand('update')
	}
	d
}

# quit
postsim <- function(x)
	NLQuit()


nl.path <- Sys.getenv('NETLOGO_PATH', NETLOGO_PATH)
model.path <- paste(getwd(), MODEL_PATH, sep='/')
parLapply(cl, 1:processors, presim, gui=FALSE, nl.path=nl.path, model.path=model.path)

result <- t(parSapply(cl, 1:NSim, simfun))

parLapply(cl, 1:processors, postsim)
stopCluster(cl)

toc <- Sys.time()
print(toc-tic)

## Process results

M <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, mean))
V <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, sd))
#MIN <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, min))
#MAX <- apply(result, 2, function(r) apply(matrix(unlist(r), nrow=NTicks), 1, max))

# standard error = standard deviation / sqrt(NSim)


library(reshape2)
X <- melt(data.frame(t=1:NTicks, M), id.vars='t')
Y1 <- melt(data.frame(t=1:NTicks, V), id.vars='t')
Z <- data.frame(X, ymin=(X$value-Y1$value), ymax=(X$value+Y1$value))

# change time scales so that the first weeks are displayed with more promeninence
Z$t[Z$t > 12] <- Z$t[Z$t > 12] + 100
Z$t[Z$t <= 12] <- (Z$t[Z$t <= 12]-1) * 100/11

output <- function(filename, width)
	pdf(paste(filename, '.pdf', sep=''),
		paper='special', width=8, height=8/width)

library(ggplot2)

output('hivsim0', 2)
p <- ggplot(Z, aes(x=t, y=value, group=variable, ymin=ymin, ymax=ymax)) +
	geom_ribbon(aes(fill=variable), alpha=0.2) +
	geom_line(aes(color=variable)) +
	xlab("") + ylab(expression(paste('# cells x ', 10^5, sep=''))) +
	scale_color_manual("Type", values=c('red','yellow','blue')) +
	scale_fill_manual("Type", values=c('red','yellow','blue')) +
	ggtitle("Cell count")
print(p)
dev.off()

